

/*
 * Class:        BrownianMotion
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */

package umontreal.iro.lecuyer.stochprocess;
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;



/**
 * This class represents a <SPAN  CLASS="textit">Brownian motion</SPAN> process 
 * <SPAN CLASS="MATH">{<I>X</I>(<I>t</I>) : <I>t</I>&nbsp;&gt;=&nbsp;0}</SPAN>,
 * sampled at times 
 * <SPAN CLASS="MATH">0 = <I>t</I><SUB>0</SUB> &lt; <I>t</I><SUB>1</SUB> &lt; <SUP> ... </SUP> &lt; <I>t</I><SUB>d</SUB></SPAN>.
 * This process obeys the stochastic differential equation
 * 
 * <P></P>
 * <DIV ALIGN="CENTER" CLASS="mathdisplay"><A NAME="eq:Brownian-motion"></A>
 * <I>dX</I>(<I>t</I>) = <I>&#956;dt</I> + <I>&#963;dB</I>(<I>t</I>),
 * </DIV><P></P>
 * with initial condition <SPAN CLASS="MATH"><I>X</I>(0) = <I>x</I><SUB>0</SUB></SPAN>,
 * where <SPAN CLASS="MATH"><I>&#956;</I></SPAN> and <SPAN CLASS="MATH"><I>&#963;</I></SPAN> are the drift and volatility parameters,
 * and 
 * <SPAN CLASS="MATH">{<I>B</I>(<I>t</I>),&nbsp;<I>t</I>&nbsp;&gt;=&nbsp;0}</SPAN> is a standard Brownian motion
 * (with drift 0 and volatility 1).
 * This process has stationary and independent increments over disjoint
 * time intervals (it is a L&#233;vy process) and the increment over an interval
 * of length <SPAN CLASS="MATH"><I>t</I></SPAN> is normally distributed with mean <SPAN CLASS="MATH"><I>&#956;t</I></SPAN> and variance 
 * <SPAN CLASS="MATH"><I>&#963;</I><SUP>2</SUP><I>t</I></SPAN>.
 * 
 * <P>
 * In this class, this process is generated using the sequential (or random walk)
 * technique:  <SPAN CLASS="MATH"><I>X</I>(0) = <I>x</I><SUB>0</SUB></SPAN> and
 * 
 * <P></P>
 * <DIV ALIGN="CENTER" CLASS="mathdisplay"><A NAME="eq:Brownian-motion-sequential"></A>
 * <I>X</I>(<I>t</I><SUB>j</SUB>) - <I>X</I>(<I>t</I><SUB>j-1</SUB>) = <I>&#956;</I>(<I>t</I><SUB>j</SUB> - <I>t</I><SUB>j-1</SUB>) + <I>&#963;</I>(t_j - t_j-1)<SUP>1/2</SUP><I>Z</I><SUB>j</SUB>
 * </DIV><P></P>
 * where 
 * <SPAN CLASS="MATH"><I>Z</I><SUB>j</SUB>&#8764;<I>N</I>(0, 1)</SPAN>.
 * 
 */
public class BrownianMotion extends StochasticProcess  {
    protected NormalGen    gen;
    protected double       mu,
                           sigma;
    // Precomputed values for standard BM
    protected double[]     mudt,
                           sigmasqrdt;



   /**
    * Constructs a new <TT>BrownianMotion</TT> with
    * parameters <SPAN CLASS="MATH"><I>&#956;</I> =</SPAN> <TT>mu</TT>, <SPAN CLASS="MATH"><I>&#963;</I> =</SPAN> <TT>sigma</TT> and initial value
    * 
    * <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>0</SUB>) =</SPAN> <TT>x0</TT>.
    * The normal variates <SPAN CLASS="MATH"><I>Z</I><SUB>j</SUB></SPAN> in will be
    * generated by inversion using <TT>stream</TT>.
    * 
    */
   public BrownianMotion (double x0, double mu, double sigma,
                          RandomStream stream) {
        this (x0, mu, sigma, new NormalGen (stream));
    }


   /**
    * Constructs a new <TT>BrownianMotion</TT> with parameters <SPAN CLASS="MATH"><I>&#956;</I> =</SPAN>
    * <TT>mu</TT>, <SPAN CLASS="MATH"><I>&#963;</I> =</SPAN> <TT>sigma</TT> and initial value 
    * <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>0</SUB>) =</SPAN> <TT>x0</TT>.
    * Here, the normal variate generator
    * {@link umontreal.iro.lecuyer.randvar.NormalGen NormalGen} is specified
    * directly instead of specifying the stream and using inversion.
    * The normal generator <TT>gen</TT> can use another method than inversion.
    * 
    */
   public BrownianMotion (double x0, double mu, double sigma, NormalGen gen) {
        this.mu    = mu;
        this.sigma = sigma;
        this.x0    = x0;
        this.gen   = gen;
    }
 

   public double nextObservation() {
        double x = path[observationIndex];
        x += mudt[observationIndex]
             + sigmasqrdt[observationIndex] * gen.nextDouble();
        observationIndex++;
        path[observationIndex] = x;
        return x;
    }


   /**
    * Generates and returns the next observation at time <SPAN CLASS="MATH"><I>t</I><SUB>j+1</SUB> =</SPAN>
    *  <TT>nextTime</TT>. It uses the previous observation time <SPAN CLASS="MATH"><I>t</I><SUB>j</SUB></SPAN> defined earlier
    * (either by this method or by <TT>setObservationTimes</TT>),
    * as well as the value of the previous observation <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>j</SUB>)</SPAN>.
    * <SPAN  CLASS="textit">Warning</SPAN>: This method will reset the observations time <SPAN CLASS="MATH"><I>t</I><SUB>j+1</SUB></SPAN>
    * for this process to <TT>nextTime</TT>. The user must make sure that
    * the <SPAN CLASS="MATH"><I>t</I><SUB>j+1</SUB></SPAN> supplied is 
    * <SPAN CLASS="MATH">&nbsp;&gt;=&nbsp;<I>t</I><SUB>j</SUB></SPAN>.
    * 
    */
   public double nextObservation (double nextTime)  {
        // This method is useful for generating variance gamma processes
        double x = path[observationIndex];
        double previousTime = t[observationIndex];
        observationIndex++;
        t[observationIndex] = nextTime;
        double dt = nextTime - previousTime;
        x += mu * dt + sigma * Math.sqrt (dt) * gen.nextDouble();
        path[observationIndex] = x;
        return x;
    }


   /**
    * Generates an observation of the process in <TT>dt</TT> time units,
    * assuming that the process has value <SPAN CLASS="MATH"><I>x</I></SPAN> at the current time.
    * Uses the process parameters specified in the constructor.
    * Note that this method does not affect the sample path of the process
    * stored internally (if any).
    * 
    */
   public double nextObservation (double x, double dt)  {
        x += mu * dt + sigma * Math.sqrt (dt) * gen.nextDouble();
        return x;
    }


   public double[] generatePath() {
        double x = x0;
        for (int j = 0; j < d; j++) {
            x += mudt[j] + sigmasqrdt[j] * gen.nextDouble();
            path[j + 1] = x;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }

   /**
    * Same as generatePath(), but a vector of uniform random numbers
    * must be provided to the method.  These uniform random numbers are used
    * to generate the path.
    * 
    */
   public double[] generatePath (double[] uniform01)  {
        double x = x0;
        for (int j = 0; j < d; j++) {
            x += mudt[j] + sigmasqrdt[j] * NormalDist.inverseF01(uniform01[j]);
            path[j + 1] = x;
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }


   public double[] generatePath (RandomStream stream) {
        gen.setStream (stream);
        return generatePath();
    }

   /**
    * Resets the parameters 
    * <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>0</SUB>) = <texttt>x0</texttt></SPAN>, 
    * <SPAN CLASS="MATH"><I>&#956;</I> = <texttt>mu</texttt></SPAN> and
    * 
    * <SPAN CLASS="MATH"><I>&#963;</I> = <texttt>sigma</texttt></SPAN> of the process.
    * <SPAN  CLASS="textit">Warning</SPAN>: This method will recompute some quantities stored internally,
    * which may be slow if called too frequently.
    * 
    */
   public void setParams (double x0, double mu, double sigma)  {
        this.x0    = x0;
        this.mu    = mu;
        if (sigma <= 0)
           throw new IllegalArgumentException ("sigma <= 0");
        this.sigma = sigma;
        if (observationTimesSet) init(); // Otherwise not needed.
    }


   /**
    * Resets the random stream of the normal generator to <TT>stream</TT>.
    * 
    */
   public void setStream (RandomStream stream)  { gen.setStream (stream); }


   /**
    * Returns the random stream of the normal generator.
    * 
    */
   public RandomStream getStream()  { return gen.getStream (); }


   /**
    * Returns the value of <SPAN CLASS="MATH"><I>&#956;</I></SPAN>.
    * 
    */
   public double getMu()  { return mu; }


   /**
    * Returns the value of <SPAN CLASS="MATH"><I>&#963;</I></SPAN>.
    * 
    */
   public double getSigma()  { return sigma; }


   /**
    * Returns the normal random variate generator used.
    * The {@link umontreal.iro.lecuyer.rng.RandomStream RandomStream}
    * used by that generator can be changed via
    * <TT>getGen().setStream(stream)</TT>, for example.
    * 
    */
   public NormalGen getGen()  { return gen; }
 

    // This is called by setObservationTimes to precompute constants
    // in order to speed up the path generation.
   protected void init() {
        super.init();
        mudt       = new double[d];
        sigmasqrdt = new double[d];
        for (int j = 0; j < d; j++) {
            double dt     = t[j+1] - t[j];
            mudt[j]       = mu * dt;
            sigmasqrdt[j] = sigma * Math.sqrt (dt);
        }
     }

} 
